Problema y Objetivo

Los sistemas para compartir bicicletas (como EcoBici en México) son una nueva generación del tradicional alquiler de bicicletas. En éstos, todo el proceso de membresía, alquiler y devolución se ha automatizado. Hoy en día, existe un gran interés en estos sistemas debido a su importante papel en la reducción del tráfico vehicular, el cuidado del medio ambiente y la reducción en los problemas de salud. Actualmente, hay alrededor de 500 programas para compartir bicicletas en todo el mundo que se componen por más de 500 mil bicicletas.

A través de estos sistemas, el usuario puede alquilar una bicicleta fácilmente desde una estación particular y regresarla en una estación que no tiene porqué coincidir con la inicial. Opuesto a otros servicios de transporte como el autobús o el metro, la duración del viaje, la posición de salida y la de llegada se registran explícitamente en estos sistemas. Esta característica convierte el sistema de intercambio de bicicletas en una red de sensores virtuales que se puede usar para detectar la movilidad en la ciudad y otras características. Por ejemplo, el proceso de alquiler de bicicletas compartidas está altamente correlacionado con las condiciones ambientales. Las condiciones climáticas, la precipitación, el día de la semana, la estación del año, la hora del día, etc. pueden afectar las conductas de alquiler.

De esta manera, se busca predecir el número de bicicletas rentadas por hora según las condiciones climáticas y ambientales y comparar los resultados obtenidos con distintos métodos vistos en clase. Para ello, se utiliza el registro histórico (años 2011 y 2012) del sistema Capital Bikeshare, Washington D.C., EE. UU.[^1] La base de datos contiene 17,379 observaciones para 17 variables; entre ellas se encuentran: la fecha y hora, estación del año, si el día es laborable o no, tipo de clima, condiciones ambientales (temperatura, humedad, etc) y número de bicicletas rentadas (variable de interés). Los métodos a utilizar para la estimación y comparación son: regresión lineal, regresión lineal con regularización, redes neuronales, bosques aleatorios y gradient boosting machine. [^1]: Los datos fueron recuperados de UCI Machine Learning Repository.

Datos y análisis descriptivo

# Se leen los datos 
datos <- read_csv("hour.csv",
    col_types =
    cols(
  instant = col_integer(),
  dteday = col_date(format = ""),
  season = col_character(),
  yr = col_integer(),
  mnth = col_character(),
  hr = col_character(),
  holiday = col_integer(),
  weekday = col_character(),
  workingday = col_integer(),
  weathersit = col_character(),
  temp = col_double(),
  atemp = col_double(),
  hum = col_double(),
  windspeed = col_double(),
  casual = col_integer(),
  registered = col_integer(),
  cnt = col_integer()))

datos_dummies <- dummy.data.frame(datos%>%data.frame, sep = "_") %>% as.tibble() %>% 
  select(-instant, -casual, -registered)

dim(datos)
## [1] 17379    17

Variables

Las variables contenidas en la base de datos son:

  • instant: índice de observaciones
  • dteday : fecha
  • season : Estación del año (1=primavera, 2=verano, 3=otoño, 4=invierno)
  • yr : año (0=2011, 1=2012)
  • mnth : mes ( 1 to 12)
  • hr : hora (0 to 23)
  • holiday : día festivo (1) o no festivo (0)
  • weekday : día de la semana (0 a 6 con 0=Domingo)
  • workingday : día laboral (1= no es fin de semana ni festivo, 0 e.o.c.)
  • weathersit : clima

    + 1 = Claro, Pocas nubes, Parcialmente nublado, Parcialmente nublado
    + 2 = Niebla + Nublado, Niebla + Nubes fragmentadas, Niebla + Pocas nubes, Niebla
    + 3 = Nieve ligera, lluvia ligera + tormenta eléctrica + nubes dispersas, lluvia ligera + nubes dispersasred clouds
    + 4 = luvia pesada + paletas de hielo + tormenta + niebla, nieve + nieblaist, Snow + Fog
  • temp : Temperatura normalizada en grados Celcius. Divididos entre 41 (máximo)
  • atemp: Sensación térmica normalizada en grados Celcius. Divididos entre 50 (máximo)
  • hum: Humedad normalizada. Divididos entre 100 (máximo)
  • windspeed: velocidad del viendo normalizada. Divididos entre 67 (máximo)
  • casual: número de usuarios casuales
  • registered: número de usuarios registrados
  • cnt: número total de bicicletas rentadas (incluye usuarios casuales y registrados)

Manipulación de datos

En primer lugar, eliminamos las variables casual y registered pues son variables que no se tendrán en la verdadera tarea de predicción. Transformarmos las variables categóricas (season, mnth, hr, weekday y weathersit ) en dummies para facilitar su manejo en algunos modelos.

# Conservamos la base con datos sin one hot encoding
datos_sin <- datos %>% 
  mutate(datetime = as.POSIXct((dteday) + hours(hr))) %>%
  select(-casual, -registered, -instant)

# Eliminamos las variables casual y registered, también creamos la variable datetime, para poder hacer una mejor visualización de los datos.
datos<-datos%>%
  select(-c(casual,registered))

# Almacenamos la fecha
fecha<-paste(datos$dteday,datos$hr,sep="/")
 
# Se asignan como factores las variables casual y registred pues  
cols<-match(c('season', 'mnth', 'hr', 'weekday','weathersit'),names(datos))
datos[cols]<-lapply(datos[,cols],as.factor)

# Se crea una dummy para cada categoría (menos una)
auxdatos<-model.matrix(~0+season+mnth+hr+weekday+weathersit,datos)

# Se eliminan las varaibles categóricas y se agregan las dummies
datos<-datos%>%
  select(-c(season, mnth, hr, weekday,weathersit))%>%
  mutate(dteday = ymd(dteday)) %>% 
  cbind(auxdatos)
head(datos) 
##   instant     dteday yr holiday workingday temp  atemp  hum windspeed cnt
## 1       1 2011-01-01  0       0          0 0.24 0.2879 0.81    0.0000  16
## 2       2 2011-01-01  0       0          0 0.22 0.2727 0.80    0.0000  40
## 3       3 2011-01-01  0       0          0 0.22 0.2727 0.80    0.0000  32
## 4       4 2011-01-01  0       0          0 0.24 0.2879 0.75    0.0000  13
## 5       5 2011-01-01  0       0          0 0.24 0.2879 0.75    0.0000   1
## 6       6 2011-01-01  0       0          0 0.24 0.2576 0.75    0.0896   1
##   season1 season2 season3 season4 mnth10 mnth11 mnth12 mnth2 mnth3 mnth4
## 1       1       0       0       0      0      0      0     0     0     0
## 2       1       0       0       0      0      0      0     0     0     0
## 3       1       0       0       0      0      0      0     0     0     0
## 4       1       0       0       0      0      0      0     0     0     0
## 5       1       0       0       0      0      0      0     0     0     0
## 6       1       0       0       0      0      0      0     0     0     0
##   mnth5 mnth6 mnth7 mnth8 mnth9 hr1 hr10 hr11 hr12 hr13 hr14 hr15 hr16
## 1     0     0     0     0     0   0    0    0    0    0    0    0    0
## 2     0     0     0     0     0   1    0    0    0    0    0    0    0
## 3     0     0     0     0     0   0    0    0    0    0    0    0    0
## 4     0     0     0     0     0   0    0    0    0    0    0    0    0
## 5     0     0     0     0     0   0    0    0    0    0    0    0    0
## 6     0     0     0     0     0   0    0    0    0    0    0    0    0
##   hr17 hr18 hr19 hr2 hr20 hr21 hr22 hr23 hr3 hr4 hr5 hr6 hr7 hr8 hr9
## 1    0    0    0   0    0    0    0    0   0   0   0   0   0   0   0
## 2    0    0    0   0    0    0    0    0   0   0   0   0   0   0   0
## 3    0    0    0   1    0    0    0    0   0   0   0   0   0   0   0
## 4    0    0    0   0    0    0    0    0   1   0   0   0   0   0   0
## 5    0    0    0   0    0    0    0    0   0   1   0   0   0   0   0
## 6    0    0    0   0    0    0    0    0   0   0   1   0   0   0   0
##   weekday1 weekday2 weekday3 weekday4 weekday5 weekday6 weathersit2
## 1        0        0        0        0        0        1           0
## 2        0        0        0        0        0        1           0
## 3        0        0        0        0        0        1           0
## 4        0        0        0        0        0        1           0
## 5        0        0        0        0        0        1           0
## 6        0        0        0        0        0        1           1
##   weathersit3 weathersit4
## 1           0           0
## 2           0           0
## 3           0           0
## 4           0           0
## 5           0           0
## 6           0           0

Separación de muestras

La separación de muestras se realizará sobre la línea de tiempo: 1. Muestra de entrenamiento: datos por hora del 1° de enero del 2011 al 30 de septiembre del 2012. 2. Muestra de validación: datos por hora del 1° de octubre del 2012 al 30 de noviembre del 2012. 3. Muestra de prueba: datos por hora del 1 de diciembre del 2012 al 31 de diciembre del 2012.

test  <- datos %>% filter(dteday >= "2012-12-01") %>% select(-instant, -dteday)
valid <- datos %>% filter(dteday >= "2012-10-01" & dteday < "2012-11-30") %>% 
  select(-instant, -dteday)
train <- datos %>% filter(dteday <  "2012-09-30") %>% select(-instant, -dteday)

test_dummies  <- datos_dummies %>% filter(dteday >= "2012-12-01") %>% select(-dteday)
valid_dummies <- datos_dummies %>% filter(dteday >= "2012-10-01" & dteday < "2012-11-30") %>% 
  select(-dteday)
train_dummies <- datos_dummies %>% filter(dteday <  "2012-09-30") %>% 
  select(-dteday)

# Datos sin one hot encoding

test_sin  <- datos_sin %>% filter(dteday >= "2012-12-01") %>% select(-datetime, -dteday)
valid_sin <- datos_sin %>% filter(dteday >= "2012-10-01" & dteday < "2012-11-30") %>%
  select(-datetime, -dteday)
train_sin <- datos_sin %>% filter(dteday <  "2012-09-30")  %>% select(-datetime, -dteday)

# Fechas almacenadas por aparte

fechas_test  <- datos_sin %>% filter(dteday >= "2012-12-01") %>% select(datetime, dteday)
fechas_valid <- datos_sin %>% filter(dteday >= "2012-10-01" & dteday < "2012-11-30") %>%
  select(datetime, dteday)
fechas_train <- datos_sin %>% filter(dteday <  "2012-09-30")  %>% select(datetime, dteday)

Una vez separadas las muestras, se estandarizan las variables numéricas (menos la respuesta) en los datos de entrenamiento y, con esas medias y varianzas, se estandarizan los datos de validación y prueba.

Visualización de datos

Primero podemos ver el número de usuarios de las bicicletas por cada hora, tomando en cuenta la temperatura y la estación del año en que se utiliza. Se puede apreciar que, como podría esperarse, hay un mayor número de usuarios entre las 7:00 y 8:00 AM, que es a la hora en que las personas se trasladan a su trabajo, y también entre 5:00 y 6:00 PM, cuando van de regreso a sus hogares. Además se puede ver que se utilizan menos las bicicletas durante las estaciones más frías, invierno y primavera.

graf <- ggplot(data = cbind(train_sin, fechas_train) )

graf +
  geom_jitter(aes(x = as.numeric(hr), y = cnt, color = temp)) +
  ggtitle("Número de usuarios de acuerdo a la hora") +
  xlab("Hora") +
  ylab("Número de usuarios") +
  scale_x_continuous(breaks=c(0:23), labels=c(0:23)) +
  facet_grid(season ~ ., 
             labeller = as_labeller(c("1" = "Primavera","2" = "Verano",
                                      "3" = "Otoño","4" = "Invierno"))) +
  theme_tufte(base_size = 14, base_family = "sans")

La estación del año podría ser de las variables importantes a la hora de que un usuario decida utilizar una bicicleta o no, ya que esperaríamos que los usuarios utilicen el servicio en las temporadas más cálidas. Esto lo podemos comprobar en la siguiente gráfica, donde se puede ver que hay un promedio de usuarios un poco más bajo durante la primavera y el invierno. Además de eso, también se aprecia un promedio menor los sábados y domingos.

ggplot() +
  geom_boxplot(data = rbind(train_sin, valid_sin), aes(x = season, y = cnt, fill =as.factor(weekday))) +
  ggtitle("Número de usuarios de acuerdo a la estación del año") +
  xlab("Estación del año") +
  ylab("Número de usuarios") +
  scale_x_discrete(labels = c("Primavera", "Verano", "Otoño", "Invierno")) +
  scale_fill_discrete(name = "Día", 
                      labels = c("Domingo", "Lunes", "Martes", 
                                 "Miércoles", "Jueves", "Viernes", "Sábado")) +
  theme_tufte(base_size = 14, base_family = "sans")

También las condiciones de lluvia y humedad pueden afectar el número de personas que deciden transportarse en bicicleta. En la siguiente gráfica se puede apreciar que esta relación sí existe para los días laborales: a mayor humedad, menor número de usuarios.

graf +
  geom_jitter(aes(x = as.numeric(hr), y = cnt, color = hum)) +
  ggtitle("Número de usuarios y humedad") + 
  xlab("Hora") +
  ylab("Número de usuarios") +
  scale_x_continuous(breaks=c(0:23), labels=c(0:23)) +
  facet_grid(workingday~.,
             labeller = as_labeller(c("1" = "Laboral","0" = "No laboral"))) +
  theme_tufte(base_size = 14, base_family = "sans")

Ajuste y Validación de modelos

En esta sección realizamos el ajuste del modelo con distintos métodos vistos en clase. Dentro de cada método, la selección de hiperparámetros se realiza con base en el menor error de validación.

Regresión Lineal

Para realizar la regresión lineal, usamos los datos sin one hot encoding.

modelo_lineal <- lm(cnt ~ . -workingday, data = train_sin)

summary(modelo_lineal)
## 
## Call:
## lm(formula = cnt ~ . - workingday, data = train_sin)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -354.14  -59.34   -6.58   49.72  437.50 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -72.296      6.794 -10.641  < 2e-16 ***
## season2       29.753      4.860   6.122 9.48e-10 ***
## season3       16.615      6.061   2.741 0.006128 ** 
## season4       42.214      6.029   7.002 2.63e-12 ***
## yr            89.030      1.767  50.393  < 2e-16 ***
## mnth10        42.722      7.958   5.369 8.05e-08 ***
## mnth11        28.400      7.808   3.637 0.000276 ***
## mnth12        29.368      6.107   4.809 1.53e-06 ***
## mnth2          5.477      3.818   1.434 0.151451    
## mnth3         22.142      4.356   5.083 3.77e-07 ***
## mnth4         21.640      6.566   3.296 0.000983 ***
## mnth5         39.269      7.095   5.535 3.17e-08 ***
## mnth6         29.884      7.430   4.022 5.79e-05 ***
## mnth7         16.840      8.508   1.979 0.047807 *  
## mnth8         36.507      8.284   4.407 1.05e-05 ***
## mnth9         61.087      7.579   8.060 8.20e-16 ***
## hr1          -17.507      5.559  -3.149 0.001639 ** 
## hr10         105.334      5.585  18.859  < 2e-16 ***
## hr11         129.542      5.626  23.026  < 2e-16 ***
## hr12         166.760      5.673  29.397  < 2e-16 ***
## hr13         162.868      5.716  28.492  < 2e-16 ***
## hr14         147.541      5.747  25.672  < 2e-16 ***
## hr15         155.252      5.761  26.949  < 2e-16 ***
## hr16         215.865      5.749  37.545  < 2e-16 ***
## hr17         370.634      5.720  64.801  < 2e-16 ***
## hr18         343.219      5.684  60.385  < 2e-16 ***
## hr19         238.309      5.627  42.349  < 2e-16 ***
## hr2          -25.709      5.581  -4.607 4.13e-06 ***
## hr20         158.771      5.594  28.383  < 2e-16 ***
## hr21         109.864      5.570  19.724  < 2e-16 ***
## hr22          72.506      5.558  13.045  < 2e-16 ***
## hr23          32.542      5.554   5.860 4.74e-09 ***
## hr3          -36.799      5.620  -6.548 6.01e-11 ***
## hr4          -39.986      5.630  -7.102 1.28e-12 ***
## hr5          -23.701      5.590  -4.240 2.25e-05 ***
## hr6           33.919      5.573   6.087 1.18e-09 ***
## hr7          164.757      5.561  29.628  < 2e-16 ***
## hr8          298.289      5.554  53.703  < 2e-16 ***
## hr9          156.913      5.561  28.217  < 2e-16 ***
## holiday      -21.005      5.292  -3.970 7.24e-05 ***
## weekday1       4.985      3.103   1.606 0.108210    
## weekday2       7.819      3.021   2.588 0.009651 ** 
## weekday3       9.348      3.024   3.091 0.001998 ** 
## weekday4       9.698      3.022   3.209 0.001335 ** 
## weekday5      13.034      3.015   4.323 1.55e-05 ***
## weekday6      12.573      2.999   4.193 2.77e-05 ***
## weathersit2   -9.653      2.031  -4.753 2.02e-06 ***
## weathersit3  -61.611      3.336 -18.469  < 2e-16 ***
## weathersit4  -62.670     57.265  -1.094 0.273804    
## temp         124.989     29.747   4.202 2.66e-05 ***
## atemp         87.224     30.725   2.839 0.004534 ** 
## hum          -82.439      5.694 -14.477  < 2e-16 ***
## windspeed    -31.182      7.290  -4.277 1.90e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 98.9 on 15134 degrees of freedom
## Multiple R-squared:  0.6915, Adjusted R-squared:  0.6904 
## F-statistic: 652.2 on 52 and 15134 DF,  p-value: < 2.2e-16

Se calcula el error cuadrático medio del modelo lineal.

preds_lineal <- predict(modelo_lineal, newdata = valid_sin)
resid_lineal <- valid_sin$cnt - preds_lineal
mse_valid_lineal <- mean(resid_lineal**2)
mse_valid_lineal
## [1] 16486.67

Se elaboran gráficas para observar el desempeño de este modelo. En la primera gráfica se puede apreciar que el modelo generalmente identifica la tendencia temporal, pero tiende a subestimar el número real de usuarios que utilizan las bicicletas. En la segunda gráfica podemos observar que los residuos del modelo se alejan mucho del 0, confirmando la subestimación. En general, se puede observar que no se logra un buen ajuste de los datos de validación.

ggplot() +
  geom_line(aes(x = fechas_valid$datetime, y = valid_sin$cnt), color = "blue", alpha = 0.5, size = 0.5) +
  geom_line(aes(x = fechas_valid$datetime, y = preds_lineal), color = "red", alpha = 0.5, size = 0.5) +
  ggtitle("Predicciones del modelo lineal y valores reales a través del tiempo") +
  ylab("Número de usuarios") +
  xlab("Fecha") +
  theme_tufte(base_size = 14, base_family = "sans")

ggplot() +
  geom_point(aes(x = fechas_valid$datetime, y = resid_lineal)) +
  geom_abline(intercept = 0, slope = 0, color = "red", size = 1.2) +
  ggtitle("Residuales del modelo lineal") +
  ylab("Residuales") +
  xlab("Fecha") +
  theme_tufte(base_size = 14, base_family = "sans")

ggplot() +
  geom_jitter(aes(x = valid$cnt, y = preds_lineal)) +
  geom_abline(intercept = 0, slope = 1, color = "red", size = 1.2) +
  labs(x = "Usuarios", y = "Usuarios predecidos") +
  ggtitle("Predicciones del modelo lineal vs. Valores reales") +
  theme_tufte(base_size = 14, base_family = "sans")

Regresión Lineal con regularización

Creamos una función que ajuste un nuevo modelo de acuerdo a cada lambda que se encuentra dentro de una lista que se utiliza como input.

modelos_lasso <- function(lambdas, alphas){
  dfs <- data.frame()
  
  x_train <- train_dummies %>% select(-cnt) %>% data.matrix()
  y_train <- train_dummies %>% select(cnt) %>% data.matrix()
  
  x_valid <- valid_dummies %>% select(-cnt) %>% data.matrix()
  y_valid <- valid_dummies %>% select(cnt) %>% data.matrix()
  
  x_test  <- test_dummies  %>% select(-cnt) %>% data.matrix()
  y_test  <- test_dummies  %>% select(cnt) %>% data.matrix()
  
  for(j in alphas){
    for(i in lambdas){
      set.seed(576327)
      mod_glm <- glmnet(x = x_train,
                        y = y_train,
                        lambda = i,
                        family = "gaussian",
                        intercept = T,
                        alpha = j)
      
      # Calculamos los residuos
      preds_train <- predict(mod_glm, newx = x_train)
      resid_train <- y_train - preds_train
      
      preds_valid <- predict(mod_glm, newx = x_valid)
      resid_valid <- y_valid - preds_valid
      
      # Calculamos el error cuadrático medio
      mse_train <- mean((resid_train)**2) %>% sqrt()
      mse_valid <- mean((resid_valid)**2) %>% sqrt()
      
      # Se crea una dataframe con info de cada modelo
      df <- data.frame(lambdas = i,
                       alphas = j,
                       MSE_Train = mse_train,
                       MSE_Valid = mse_valid
                       )
      # Agregamos el dataframe anterior a dfs
      dfs <- rbind(dfs, df)
      
    }
    
  }
  # Tomamos la lambda que minimiza el error cuadrático medio de valid
  lambda_min <- dfs$lambdas[which.min(dfs$MSE_Valid)]
  alpha_min <- dfs$alphas[which.min(dfs$MSE_Valid)]
  
  # Corremos el modelo con los hiperparámetros anteriores
  glm_mejor <- glmnet(x = x_train,
                        y = y_train,
                        lambda = lambda_min,
                        family = "gaussian",
                        intercept = T,
                        standardize = F,
                        alpha = alpha_min)
  
  # Generamos los residuos de validación del mejor modelo
  preds <- predict(glm_mejor, newx = x_valid)
  resids <- y_valid - preds
  
  return(list(modelos = dfs, 
              mejor = glm_mejor, 
              x_train = x_train,
              x_valid = x_valid,
              x_test  = x_test,
              y_train = y_train,
              y_valid = y_valid,
              y_test  = y_test,
              preds = preds,
              resids = resids))
}

Corremos 222 modelos, tanto LASSO como ridge, y los ordenamos por error cuadrático medio de validación.

modelos_lass <- modelos_lasso(lambdas = exp(seq(-10,10,.05)), alphas = c(0,1))
modelos_lass[[1]] %>% arrange(MSE_Valid) %>% head(10) %>% knitr::kable()
lambdas alphas MSE_Train MSE_Valid
4.54e-05 0 98.73236 128.3908
4.77e-05 0 98.73236 128.3908
5.02e-05 0 98.73236 128.3908
5.27e-05 0 98.73236 128.3908
5.55e-05 0 98.73236 128.3908
5.83e-05 0 98.73236 128.3908
6.13e-05 0 98.73236 128.3908
6.44e-05 0 98.73236 128.3908
6.77e-05 0 98.73236 128.3908
7.12e-05 0 98.73236 128.3908

Graficamos las lambdas para cualesquier alpha. La línea roja vertical indica la lambda con la que se minimiza el error cuadrático medio de validación.

ggplot() +
  facet_grid(alphas~ .) +
  geom_point(data = modelos_lass[[1]], aes(x = log(lambdas), y = (MSE_Valid), color =alphas), size=0.7) +
  geom_vline(xintercept = modelos_lass[[1]]$lambda[which.min(modelos_lass[[1]]$MSE_Valid)] %>% log(),
             color = "red") +
  ylab("MSE Validación") +
  theme_tufte(base_size = 14, base_family = "sans")

La lambda que minimiza el error cuadrático medio de validación es muy pequeña, por lo que los resultados de este modelo se parecen mucho a los de la regresión lineal de la sección anterior.

ggplot() +
  geom_line(aes(x = fechas_valid$datetime, y = modelos_lass$y_valid), color = "blue") +
  geom_line(aes(x = fechas_valid$datetime, y = modelos_lass$preds), color = "red") +
  ylab("Número de usuarios") +
  xlab("Fecha") +
  theme_tufte(base_size = 14, base_family = "sans")

ggplot() +
  geom_point(aes(x = fechas_valid$datetime, y = modelos_lass$resids)) +
  geom_abline(intercept = 0, slope = 0, color = "red", size = 1.2) +
  ggtitle("Residuales del modelo LASSO") +
  ylab("Residuales") +
  xlab("Fecha") +
  theme_tufte(base_size = 14, base_family = "sans")

ggplot() +
  geom_point(aes(x = valid$cnt, y = modelos_lass$preds)) +
  geom_abline(intercept = 0, slope = 1, color = "red", size = 1.2) +
  labs(x = "Usuarios", y = "Usuarios predecidos") +
  ggtitle("Predicciones de LASSO vs. Valores reales") +
  theme_tufte(base_size = 14, base_family = "sans")

Redes Neuronales

Ahora ajustamos un modelo de redes neuronales.

modelos_redes <- function(n.units, dropout_rate = 0){
  x_train <- modelos_lass$x_train
  y_train <- modelos_lass$y_train
  x_valid <- modelos_lass$x_valid
  y_valid <- modelos_lass$y_valid
  
  models <- list()
  history <- list()
  
  for(j in dropout_rate){
    for(i in n.units){
      set.seed(576327)
      model <- keras_model_sequential()
      
      model %>% 
        layer_dense(units = i, 
                    activation = 'relu',
                    kernel_initializer=initializer_random_uniform(minval=-0.5,maxval=0.5),
                    input_shape = c(ncol(x_train))) %>% 
        layer_dropout(rate = j) %>% 
        layer_dense(units = floor(i/2),
                    activation = 'relu',
                    kernel_initializer=initializer_random_uniform(minval=-0.5,maxval=0.5)) %>%
        layer_dropout(rate = j) %>%
        layer_dense(units = 1)
      
      model %>% compile(
        loss = 'mse',
        optimizer = "adam",
        metrics = c('mse','mae')
        )
      
      iteraciones <- model %>% fit(
        x_train, 
        y_train, 
        epochs = 1000, 
        batch_size= nrow(x_train),
        verbose = 0,
        validation_data = list(x_valid,y_valid)
        )
      
      models <- c(models, model)
      history <- c(history, iteraciones)
    }
  }
  return(list(models = models,
              histories = history))
    
  
}

Calculamos los residuos de la red neuronal.

modelos_neuronales <- modelos_redes(c(124, 256, 512), c(0, 0.3, 0.5))
preds_neuronales <- modelos_neuronales$models[[9]] %>% predict_on_batch(modelos_lass$x_valid)
resid_neuronales <- modelos_lass$y_valid - preds_neuronales

Graficamos los resultados anteriores. La red neuronal parece dar mucho mejores resultados que los modelos intentados previamente, aunque en la segunda gráfica se puede apreciar que este modelo no parece generalizar correctamente algunas observaciones entre el 15 y el 30 de noviembre.

ggplot() +
  geom_line(aes(x = fechas_valid$datetime, y = modelos_lass$y_valid), color = "blue", alpha = 0.5, 
            size = 0.5) +
  geom_line(aes(x = fechas_valid$datetime, y = preds_neuronales), color = "red", alpha = 0.5, size = 0.5) +
  ggtitle("Predicciones de red neuronal y valores reales a través del tiempo") +
  ylab("Número de usuarios") +
  xlab("Fecha") +
  theme_tufte(base_size = 14, base_family = "sans")

ggplot() +
  geom_point(aes(x = fechas_valid$datetime, y = resid_neuronales)) +
  geom_abline(intercept = 0, slope = 0, color = "red", size = 1.2) +
  ggtitle("Residuales de red neuronal") +
  ylab("Residuales") +
  xlab("Fecha") +
  theme_tufte(base_size = 14, base_family = "sans")

ggplot() +
  geom_point(aes(x = valid$cnt, y = preds_neuronales)) +
  geom_abline(intercept = 0, slope = 1, color = "red", size = 1.2) +
  labs(x = "Usuarios", y = "Usuarios predecidos") +
  ggtitle("Predicciones de red neuronal vs. Valores reales") +
  theme_tufte(base_size = 14, base_family = "sans")

Bosque Aleatorio

Construimos una función para correr modelos con diferentes parámetros:

bosques <- function(train, valid, ntrees){
  mtry <- c(1:12)
  forest <- list()
  dfs <- data.frame()
  for(i in mtry){
    set.seed(576327)
    bosque <- randomForest(cnt ~ ., 
                           data = train_sin, 
                           ntree = ntrees, 
                           mtry = i, 
                           importance=TRUE)
    
    # Encontramos las predicciones que arroja el modelo
    # probs_train <- predict(bosque, newdata = train)
    probs_valid <- predict(bosque, newdata = valid_sin)
    
    # Calculamos el error cuadrático medio
    mse_valid <- ((valid_sin$cnt - probs_valid)^2) %>% mean()
    
    df <- data.frame(mtry = i, 
                     ntrees = ntrees, 
                     MSE_Valid = mse_valid 
                     )
    
    dfs <- rbind(dfs, df)
      
    forest <- c(forest, list(bosque))
  }
  # return(forest)
  return(dfs)
}
number_trees <- function(ntrees = c(100, 500)){
  dfs <- data.frame()
  for(i in ntrees){
    forest <- bosques(train, valid, i)
    
    dfs <- rbind(dfs, forest)
  }
  return(dfs)
}

Corremos los modelos:

modelos_forests <- number_trees()

Los ordenamos de acuerdo al error cuadrático medio de validación.

modelos_forests %>% arrange(MSE_Valid) %>% head(10) %>% knitr::kable()
mtry ntrees MSE_Valid
10 500 5116.339
9 500 5174.625
9 100 5189.864
11 500 5193.212
8 100 5199.832
12 500 5208.819
12 100 5232.009
8 500 5259.124
11 100 5290.158
10 100 5312.424

Corremos un nuevo modelo utilizando los datos del mejor de los bosques aleatorios obtenido.

indice_forests <- which.min(modelos_forests$MSE_Valid)

set.seed(576327)
bosque_mejor <- randomForest(cnt ~ ., 
                             data = train_sin, 
                             ntree = modelos_forests[indice_forests,2], 
                             mtry = modelos_forests[indice_forests,1], 
                             importance=TRUE)
bosque_mejor
## 
## Call:
##  randomForest(formula = cnt ~ ., data = train_sin, ntree = modelos_forests[indice_forests,      2], mtry = modelos_forests[indice_forests, 1], importance = TRUE) 
##                Type of random forest: regression
##                      Number of trees: 500
## No. of variables tried at each split: 10
## 
##           Mean of squared residuals: 1467.839
##                     % Var explained: 95.35

Habiendo obtenido el mejor modelo de bosques aleatorios, ahora podemos ver cómo se desempeña. Para lograr lo anterior, se elaboran tres gráficas que aparecen a continuación. La primera de ella parece mostrarnos que se logró un buen ajuste para los datos de valdiación. En la segunda podemos observar los residuales y podemos percatarnos de que entre el 15 y el 30 de noviembre hay un periodo de tiempo en el que el modelo subestima la cantidad de usuarios que utilizan las bicicletas. Finalmente, en la tercera gráfica podemos observar que el modelo tiende a subestimar el número de usuarios cuando el número real de estos es mayor a 500.

preds_bosque <- predict(bosque_mejor, newdata = valid_sin)

ggplot() +
  geom_line(aes(x = fechas_valid$datetime, y = valid_sin$cnt), color = "blue", alpha = 0.5, size = 0.5) +
  geom_line(aes(x = fechas_valid$datetime, y = preds_bosque), color = "red", alpha = 0.5, size = 0.5) +
  ggtitle("Predicciones de bosques aleatorios y valores reales a través del tiempo") +
  ylab("Número de usuarios") +
  xlab("Fecha") +
  theme_tufte(base_size = 14, base_family = "sans")

ggplot() +
  geom_point(aes(x = fechas_valid$datetime, y = valid$cnt - preds_bosque)) +
  geom_abline(intercept = 0, slope = 0, color = "red", size = 1.2) +
  ylab("Residuales") +
  xlab("Fecha") +
  ggtitle("Residuales de bosques aleatorios") +
  theme_tufte(base_size = 14, base_family = "sans")

ggplot() +
  geom_point(aes(x = valid$cnt, y = preds_bosque)) +
  geom_abline(intercept = 0, slope = 1, color = "red", size = 1.2) +
  labs(x = "Usuarios", y = "Usuarios predecidos") +
  ggtitle("Predicciones de bosques aleatorios vs. Valores reales") +
  theme_tufte(base_size = 14, base_family = "sans")

Gradient Boosting Machine

Se modela el problema con GBM probando distintos valores para los hiperparámetros n.trees, interaction.depth, shrinkage y bag.fraction. En particular se consideran los siguientes valores de los hiperparámetros:

  • n.trees: 100, 250, 500, 750, 1000
  • interaction.depth: 1, 2, 3, 4
  • shrinkage: 0.001,0.005, 0.01, 0.05, 0.1, 0.5, 1
  • bag.fraction: 0.25, 0.5, 0.75, 1

Dado que gbm utiliza las primeras train.fraction observaciones para entrenamiento y las siguientes 1-train.fraction observaciones para hacer la validación, debemos ingresar los datos como una matriz de la forma \(data=[train, valid]'\) y especificar el porcentaje que ocupa la muestra de entrenamiento.

#Función gbm
ajustar_boosting<-function(data.frame){
  
  out_fun<-function(n.trees=500, interaction.depth=1,
                    shrinkage=0.1, bag.fraction=1,...){
    
    set.seed(576327)
    mod_boosting<-gbm(cnt~.,
                    data=data.frame,
                    distribution='gaussian', # Reg. logística
                    n.trees=n.trees, #no. arboles
                    interaction.depth =interaction.depth, #Grado de interacción variables
                    shrinkage=shrinkage, # Tasa de aprendizaje
                    bag.fraction=bag.fraction, # % de datos de entrenamiento
                    train.fraction = nrow(train_sin)/nrow(data.frame), # Procentaje para entrenamiento
                    keep.data=TRUE)
    mod_boosting
    
  }
  
  out_fun
  
}


# Evaluación Modelos 
Eval_gbm<-function(data){
  
  fun_eval<-function(gbm_ajust){
    
    # se seleccionan las variables indpe. 
    data_x<-data%>% select(-cnt)
    
    # Se obtienen las prediciones
    cnt_hat<-predict.gbm(gbm_ajust,newdata=data_x,n.trees=gbm_ajust$n.trees,type='response')
    
    # Se obtiene el ECM
    sqrt(mean((data$cnt-cnt_hat)^2))
  }
}

En total se corren 560 modelos (correspondientes a todas las posibles coibnaciones de hiperparámetros) y se elige el modelo que tenga menor error de validación. La siguiente tabla muestra los resultados obtenidos.

# Se obtienen todas las combinaciones posibles de hiperparámetros
# params<-list(n.trees=c(100,250,500,750,1000),
#                   interaction.depth=c(1,2,3,4),
#                   shrinkage=c(0.001,0.005,0.01,0.05,0.1,0.5,1),
#                   bag.fraction=c(0.25,0.5,0.75,1))%>%expand.grid

params<-list(n.trees=c(100,250,500,1000),
                  interaction.depth=c(1,2,3),
                  shrinkage=c(0.001,0.01,0.1,1),
                  bag.fraction=c(0.25,0.5,0.75,1))%>%expand.grid

# params<-list(n.trees=c(100,250),
#                   interaction.depth=c(1),
#                   shrinkage=c(0.1),
#                   bag.fraction=c(0.25))%>%expand.grid


# Se juntan las muestras de entrenamiento y validación. 
data_train_gbm<-rbind(train,valid)


# %%%%%%%%%%%%%%%%%%%%%% Calibración hiperparámetros modelo cocina %%%%%%%%%%%%%%%%%%%%%%

# Se evalúan la funciones ajustar_boosting en los datos de entrenamiento y validacion
gbm_hiper<-ajustar_boosting(data_train_gbm)

# Se evalua la función Eval_gbm en la muestra de entrenamiento y de validación.
Eval_train_gbm<-Eval_gbm(train)
Eval_valid_gbm<-Eval_gbm(valid)

# Se corre un modelo gbm para cada combinación de hiperparámetros.
modelos_gbm<-params%>%
  mutate(modelo=pmap(.,gbm_hiper))%>%
  mutate(Err_Train=unlist(lapply(modelo,Eval_train_gbm)))%>%
  mutate(Err_Valid=unlist(lapply(modelo,Eval_valid_gbm)))%>%
  select(n.trees,interaction.depth,shrinkage,bag.fraction,Err_Train,Err_Valid)%>%
  arrange(Err_Valid)

# Se imprime la tabla de modelos. 
modelos_gbm %>% head(10) %>% knitr::kable(digits = 2)
n.trees interaction.depth shrinkage bag.fraction Err_Train Err_Valid
1000 3 1 1.00 31.44 68.66
500 3 1 1.00 37.06 70.36
500 3 1 0.75 38.38 71.53
250 3 1 0.75 43.93 72.15
1000 2 1 1.00 43.48 72.39
1000 3 1 0.75 34.00 72.45
1000 2 1 0.75 44.77 74.26
500 2 1 1.00 46.91 74.28
250 3 1 1.00 43.87 74.91
1000 2 1 0.50 48.31 75.05

Se toman los hiperparámetros de el “mejor modelo” bajo el método gbm y se obtienen los siguientes resultados.

# Se corre el modelo gbm con los hiperparámetros que dan el mejor modelo
mod_gbm_best<-gbm_hiper(n.trees=modelos_gbm$n.trees[1], 
                     interaction.depth=modelos_gbm$interaction.depth[1],
                     shrinkage=modelos_gbm$shrinkage[1],
                     bag.fraction=modelos_gbm$bag.fraction[1])

# Se obtiene el desempeño del modelo en cada árbol agregado. 
graf_eval<-data.frame(train=sqrt(mod_gbm_best$train.error),
                      valid=sqrt(mod_gbm_best$valid.error),
                      n_arbol=1:length(mod_gbm_best$train.error))%>%
  gather(tipo,error,-n_arbol)

# Se grafica el desemeño del modelo. 
ggplot(graf_eval,aes(x=n_arbol,y=error,colour=tipo,group=tipo))+
  geom_line() +
  theme_tufte()

Se obtiene el ajuste del mejor modelo.

# Se obtienen las predicciones de entrenamiento y prueba
cnt_Train_gbm<-predict.gbm(mod_gbm_best,newdata=train%>%select(-cnt),n.trees=mod_gbm_best$n.trees,type='response')
cnt_Valid_gbm<-predict.gbm(mod_gbm_best,newdata=valid%>%select(-cnt),n.trees=mod_gbm_best$n.trees,type='response')

# Se grafican los datos .
data_ajuste_gbm<-data.frame(Muestra=c(rep("Train",nrow(train)),rep("Valid",nrow(valid))),
                                 Fecha=fechas_train %>% rbind(fechas_valid),
                                 cnt_Obs=c(train$cnt,valid$cnt),
                                 cnt_Ajust=c(cnt_Train_gbm,cnt_Valid_gbm))

ggplot() +
  geom_line(aes(x = fechas_valid$datetime, y = valid_sin$cnt), color = "blue", alpha = 0.5, size = 0.5) +
  geom_line(aes(x = fechas_valid$datetime, y = cnt_Valid_gbm), color = "red", alpha = 0.5, size = 0.5) +
  ggtitle("Predicciones de bosques aleatorios y valores reales a través del tiempo") +
  theme_tufte(base_size = 14, base_family = "sans")

ggplot() +
  geom_point(aes(x = fechas_valid$datetime, y = valid_sin$cnt - cnt_Valid_gbm)) +
  geom_abline(intercept = 0, slope = 0, color = "red", size = 1.2) +
  ggtitle("Residuales de bosques aleatorios") +
  theme_tufte(base_size = 14, base_family = "sans")

ggplot(data_ajuste_gbm,aes(x=cnt_Obs,y=cnt_Ajust,colour=Muestra)) +
  geom_point() +
  geom_abline(slope=1,intercept=0,lwd=1.5) +
  theme_tufte()

Comparación de modelos

En esta sección se lleva a cabo la comparación final de los modelos utilizando la muestra de prueba. El modelo que posea el menor error cuadrático medio de prueba será el que se seleccionará como el mejor modelo.

Regresión Lineal

Calculamos el error cuadrático medio del modelo lineal.

preds_prueba_lineal <- predict(modelo_lineal, newdata = test_sin)
resids_prueba_lineal <- test_sin$cnt - preds_prueba_lineal
MSEp_lineal <- (resids_prueba_lineal**2) %>% mean()
MSEp_lineal
## [1] 11706.01

Ahora podemos ver cómo ajusta el modelo a los datos de prueba.

ggplot() +
  ggtitle("Ajuste del modelo lineal") +
  ylab("Número de usuarios") +
  xlab("Fecha") +
  geom_line(aes(x = fechas_test$datetime, y = modelos_lass$y_test, colour = "Observado")) +
  geom_line(aes(x = fechas_test$datetime, y = preds_prueba_lineal, colour = "Ajustado")) +
  scale_colour_manual("", 
                      breaks = c("Observado", "Ajustado"),
                      values = c("blue", "red")) +
  theme_tufte(base_size = 14, base_family = "sans")

Regresión Lineal con regularización

Calculamos el error cuadrático medio del modelo lineal con regularización. En este caso la regularización afectó de manera negativa al modelo.

preds_prueba_lasso <- predict(modelos_lass$mejor, newx = modelos_lass$x_test)
resids_prueba_lasso <- modelos_lass$y_test - preds_prueba_lasso
MSEp_lasso <- (resids_prueba_lasso**2) %>% mean()
MSEp_lasso
## [1] 11708.74

Ahora vemos cómo ajusta el modelo a los datos de prueba.

ggplot() +
  ggtitle("Ajuste del modelo LASSO") +
  ylab("Número de usuarios") +
  xlab("Fecha") +
  geom_line(aes(x = fechas_test$datetime, y = modelos_lass$y_test, colour = "Observado")) +
  geom_line(aes(x = fechas_test$datetime, y = preds_prueba_lasso, colour = "Ajustado")) +
  scale_colour_manual("", 
                      breaks = c("Observado", "Ajustado"),
                      values = c("blue", "red")) +
  theme_tufte(base_size = 14, base_family = "sans")

Red Neuronal

Calculamos el error cuadrático medio del modelo de redes neuronales.

preds_prueba_neuronales <- modelos_neuronales$models[[9]] %>% predict_on_batch(modelos_lass$x_test)
resids_prueba_neuronales <- modelos_lass$y_test - preds_prueba_neuronales
MSEp_neuronales <- (resids_prueba_neuronales**2) %>% mean()
MSEp_neuronales
## [1] 4569.598

Ahora observamos cómo ajusta el modelo a los datos de prueba.

ggplot() +
  ggtitle("Ajuste de la red neuronal") +
  ylab("Número de usuarios") +
  xlab("Fecha") +
  geom_line(aes(x = fechas_test$datetime, y = modelos_lass$y_test, colour = "Observado")) +
  geom_line(aes(x = fechas_test$datetime, y = preds_prueba_neuronales, colour = "Ajustado")) +
  scale_colour_manual("", 
                      breaks = c("Observado", "Ajustado"),
                      values = c("blue", "red")) +
  theme_tufte(base_size = 14, base_family = "sans")

Bosques Aleatorios

Calculamos el error cuadrático medio del modelo de bosques aleatorios.

preds_prueba_bosques <- predict(bosque_mejor, newdata = test_sin)
resids_prueba_bosques <- test_sin$cnt - preds_prueba_bosques
MSEp_bosques <- resids_prueba_bosques**2 %>% mean()
MSEp_bosques
## [1] 5190.826

Ahora podemos ver cómo ajusta el modelo a los datos de prueba.

ggplot() +
  ggtitle("Ajuste de bosques aleatorios") +
  ylab("Número de usuarios") +
  xlab("Fecha") +
  geom_line(aes(x = fechas_test$datetime, y = modelos_lass$y_test, colour = "Observado")) +
  geom_line(aes(x = fechas_test$datetime, y = preds_prueba_bosques, colour = "Ajustado")) +
  scale_colour_manual("", 
                      breaks = c("Observado", "Ajustado"),
                      values = c("blue", "red")) +
  theme_tufte(base_size = 14, base_family = "sans")

Gradient Boosting Machine

Se obtiene el ajuste del mejor modelo gbm para los datos de prueba.

# Se evalua la función Eval_gbm en la muestra de prueba
Eval_test_gbm<-Eval_gbm(test)


# Se obtienen las predicciones de cnt para la muestra de prueba
cnt_Test_gbm<-predict.gbm(mod_gbm_best,newdata=test %>% select(-cnt),n.trees=mod_gbm_best$n.trees, type='response')

# Se obtiene el error de prueba
TestErr_gbm<-Eval_test_gbm(mod_gbm_best)
TestErr_gbm**2
## [1] 4658.615

Ahora observamos cómo ajusta el modelo a los datos de prueba.

# Se grafican las observaciones
data_Pred_gbm<-data.frame(Fecha=fechas_test,
                            cnt_Obs=test$cnt,
                            cnt_Ajust=cnt_Test_gbm)

ggplot(data_Pred_gbm,aes(x=cnt_Obs,y=cnt_Ajust))+geom_point(col='olivedrab')+
  geom_abline(slope=1,intercept=0,lwd=1.5) +
  theme_tufte(base_size = 14, base_family = "sans")

ggplot() +
  ggtitle("Ajuste de Gradient Boosting Machine") +
  ylab("Número de usuarios") +
  xlab("Fecha") +
  geom_line(aes(x = fechas_test$datetime, y = modelos_lass$y_test, colour = "Observado")) +
  geom_line(aes(x = fechas_test$datetime, y = cnt_Test_gbm, colour = "Ajustado")) +
  scale_colour_manual("", 
                      breaks = c("Observado", "Ajustado"),
                      values = c("blue", "red")) +
  theme_tufte(base_size = 14, base_family = "sans")

Conclusión

En la siguiente se pueden observar los errores cuadráticos medios de prueba de los modelos anteriores. La red neuronal es el modelo con el menor error cuadrático medio de prueba, por lo tanto es el modelo que ha sido seleccionado por el equipo. Además, en la sección anterior se puede ver que las predicciones de este modelo son las que mejor se ajustan a los datos observados. A pesar de esto, el modelo falla al predecir el número de usuarios en los últimos días de diciembre. Esto se puede deber a que en ese periodo la gente suele utilizar menos la bicicleta por encontrarse en vacaciones navideñas. En estudios posteriores que cuenten con información de más años se podría incorporar esta variable, independientemende de la variable holiday, para verificar que esa es la razón por la que el modelo no predice correctamente los valores de esas fechas.

data.frame(Modelo = c("Lineal", "LASSO", "Red Neuronal", "Bosque Aleatorio", "GBM"),
           `MSE prueba` = c(MSEp_lineal, MSEp_lasso, MSEp_neuronales, 
                            MSEp_bosques, TestErr_gbm**2)) %>% 
  knitr::kable(digits = 0)
Modelo MSE.prueba
Lineal 11706
LASSO 11709
Red Neuronal 4570
Bosque Aleatorio 5191
GBM 4659